Hospital Entry Simulation

Background:

Let’s take a simple model of a three department hospital, with an Emergency Department (ED), an Intensive Care Unit (ICU), and a department for General (Gen) Admission for less severe injuries. We can perform as simulation over the course of a year and see how many patients we can cure, while minimizing the amount we have to transfer, and who pass away.

Simulations are generally function-heavy. This is because it allows us perform the same operation without unnecessary copy-pasting of code. So let’s put our initial set-up of variables, which have been more-or-less arbitrarily chosen, into a function

library(dplyr)
library(ggplot2)
library(gganimate)
initial_setup   <<- function(){
  ed_beds.i       <<- 35 # Emergency Department beds available
  full_ed_beds.i  <<- 0  # Occupied Emergency Department beds
    
  icu_beds.i      <<- 50 # ICU Department beds available
  full_icu_beds.i <<- 0   # Occupied ICU Department beds
  
  gen_beds.i      <<- 50  # General Department beds available
  full_gen_beds.i <<- 0   # Occupied General Department beds
  
  total_patients  <<- 5000
  final_day       <<- 365 # final day of our simulation
  day.v           <<- seq(1:final_day)
}
initial_setup()

We’ll need a data frame to hold all of our patient data. For the sake of efficiency, we generate all of this data up front.

patient.df_setup <- function(){
  
  entry_day.v  <- sort(sample(0:final_day, total_patients, replace=TRUE)) # The day that the patient is injured
  severity.v   <- runif(n = total_patients, min = .1, max = 1) # The level of severity of patient injury
  location.v   <- rep("Home", total_patients)    # Current location of patients
  status.v     <- rep("Healthy", total_patients) # other catagories will be Injured, Released, Deceased, Cured
  placement.v  <- ifelse(severity.v >= 0.8,  "ED",
                  ifelse(severity.v > 0.6, "ICU",
                                           "Gen")) # The initial placement of patient based on severity
  
  df <- data.frame("Entry_day" = entry_day.v, "Severity" = severity.v, 
                   "Placement" = placement.v, "Location" = location.v,
                   "Status"    = status.v)
  df <- df %>% mutate(Placement = as.character(Placement), 
                                  Location  = as.character(Location),
                                  Status    = as.character(Status))
     
  df <- df %>% arrange(Entry_day)
  return(df)                         
}
patient.df <- patient.df_setup()

Now before we actually begin our simulation, let’s break it into pieces for the sake of maintainable code. This may look a bit strange at first, but it will all come together as we move toward the end.

The first function we’ll look at handles new patients entering the hospital. They are sorted into their appropriate department, and their status and location are adjusted appropriately.

entrance <- function(df, day){ # Entrance of new patients for the day
  
  df <- df %>%  mutate(Placement = ifelse(Severity >= 0.8,  "ED", ifelse(Severity > 0.6, "ICU", "Gen")),
                       Status    = ifelse((Entry_day == day) & (Location == "Home"),"Injured", Status),
                       Location  = ifelse((Status == "Injured") & (Location == "Home"),Placement, Location)) %>%
                arrange(Entry_day,Severity)
  return(df)
  }

Our second function is going to adjust our patient’s “Severity” value on a daily basis. The patient’s condition can improve, stay about the same, or deteriorate. The beta distribution has the appropriate qualities to fit our model, but it’s important to note that you could easily defend using many other distributions for this.

When we call this function later, we’re going to stay optimistic and say that patients will generally improve during their stay in the hospital.

improvements <- function(df, skew = "None"){ # Adjusts patient severity by a random value drawn from beta distribution
  # Beta distribution parameters
  if      (skew == "Left" ){param1 = 40 # Left skewed --- Patient will generally improve
                            param2 = 50} 
  else if (skew == "Right"){param1 = 50 # Right skewed --- Patient will generally deteriorate
                            param2 = 40}
  else                     {param1 = 45 # No skew --- Patient has an equal probability to improve/deteriorate
                            param2 = 45}
                  
  df <- df %>% mutate(Severity = ifelse( Location == "ED",  Severity + (rbeta(1, param1, param2) - .5) * 1.5,
                                 ifelse( Location == "ICU", Severity + (rbeta(1, param1, param2) - .5) * 1.25,
                                 ifelse( Location == "Gen", Severity + (rbeta(1, param1, param2) - .5), 
                                         Severity))),
                      Location = ifelse( Severity <= 0, "Home",  # Patient is cured and sent home.
                                 ifelse( Severity >= 1, "Morgue",# Patient passes away and is sent to the morgue
                                         Location)),
                      Status   = ifelse( Severity <= 0, "Cured", 
                                 ifelse( Severity >= 1, "Deceased", 
                                         Status)))
return(df)
}

The third and fourth functions are simply going to see how many patients are in each location, and which status they are presenting. At the end of each day, these are put into a data frame for later visualization.

tally_beds <- function(df){ # gets current number of full beds
  
  old_deceased.i    <<- deceased.i
  old_transferred.i <<- transferred.i
  
  full_ed_beds.i  <<- tally(df %>% filter(Location == "ED")) %>% as.integer()
  full_icu_beds.i <<- tally(df %>% filter(Location == "ICU")) %>% as.integer()
  full_gen_beds.i <<- tally(df %>% filter(Location == "Gen")) %>% as.integer()
  deceased.i      <<- tally(df %>% filter(Status == "Deceased")) %>% as.integer()
  cured.i         <<- tally(df %>% filter(Status == "Cured")) %>% as.integer()
  healthy.i       <<- tally(df %>% filter(Status == "Healthy")) %>% as.integer()
  transferred.i   <<- tally(df %>% filter(Status == "Transferred")) %>% as.integer()
  healhty.i       <<- tally(df %>% filter(Status == "Healthy")) %>% as.integer()
}

day_end_tally <- function(df, day){
  tally_beds(df)
  
  full_beds.df <<- rbind(full_beds.df, data.frame(Day         = day,
                                                  ED          = full_ed_beds.i,
                                                  ICU         = full_icu_beds.i, 
                                                  Gen         = full_gen_beds.i,
                                                  Healthy     = healthy.i,
                                                  Cured       = cured.i,
                                                  Deceased    = deceased.i,
                                                  Transferred = transferred.i))
}

The rerank function simply gives our lowest severity patients (per group) a rank closer to zero. Why do we need that? Because the hospital, when faced with space constraints, will try to release or transfer their lowest severity patients (with some caveats below).

rerank <- function (df){ # Ranks patient Severity by location --- lowest Severity gets lowest number rank
  df <- df %>% group_by(Location) %>%
               mutate(Rank = ifelse(Location %in% c("ED","ICU","Gen"), rank(Severity, ties.method = "first"),NA)) %>%
               ungroup()
  tally_beds(df)
  return(df)
  }

Below is our main function. It simulates one day by having new patients enter, transferring patients as needed, tallying up the beds, adjusting our patient’s severity each day, and then the function will recursively call itself until it hits the end of the simulation.

For the sake of the model, we only allow transfers to other departments if their severity is on the border between the established points. For example, a patient must be between severity level 0.8 and 0.85 to be considered for transfer from the Emergency Department to the ICU. If the hospital still can’t handle the number of patients, they are transferred to another hospital.

next_day <- function(df, day, final_day){  # Recursive function --- reassigns patients to new area if beds full
  
  df <- entrance(df, day) 
  df <- rerank(df)
  
  # reassign ED patients to ICU department if ED is full if possible
  # if patient is too severe for ICU - they are transfered to another hospital
  df <- df %>% mutate(Location = ifelse((Location == "ED")  & (Rank <= (full_ed_beds.i - ed_beds.i)),
                                 ifelse(Severity < 0.85,"ICU","Another Hospital"), Location))  
  df <- rerank(df)
  
  
  # reassign ICU patients to Gen departmentt if ICU is full
  # if patient is too severe for General admission - they are transfered to another hospital
  df <- df %>% mutate(Location = ifelse((Location == "ICU") & (Rank <= (full_icu_beds.i - icu_beds.i)),
                                 ifelse(Severity < 0.65,"Gen","Another Hospital"), Location))
  
  df <- rerank(df)
  
  # relases patients if beds are full and their severity is low enough
  # otherwise transfers them to another hospital
  df <- df %>% mutate(Location = ifelse((Location == "Gen") & (Rank <= (full_gen_beds.i - gen_beds.i)),
                                 ifelse(Severity < 0.20,"Released","Another Hospital"), Location),
                      Status   = ifelse(Location == "Another Hospital", "Transferred", Status),
                      Status   = ifelse(Location == "Released", "Released", Status))
  day_end_tally(df,day)
  df <- improvements(df, skew = "Left") # Adjusts patients severity for the day
  
  day = day + 1
  if(day <= final_day){  
    next_day (df, day, final_day) # Move on to next day if not past the final day
    }
  else {return(df)}
}

And our last function doing any analytic work will call the functions above and give us our data in a nice format for our simulation visualization.

bed_count <- function(){
  full_beds.df <<- data.frame(Day         = integer(),
                              ED          = integer(),
                              ICU         = integer(), 
                              Gen         = integer(),
                              Healthy     = integer(),
                              Cured       = integer(),
                              Deceased    = integer(),
                              Transffered = integer(),
                              Healthy     = integer())
  
  deceased.i    <<- 0
  transferred.i <<- 0
  
  patient.df <<- next_day(patient.df, 0, final_day)
  tally_beds(patient.df)
}
bed_count()
full_beds.df

Above, in tabular form, we can see how our departments and statuses of our patients change each day. But there is only so much data that can be absorbed through a table. In order for these data to be meaningful, we have to understand it on a deeper level than can be conveyed by a static number on a page.

And for that, we use visualizations. Using the gganimate package, and the data generated above, we can show how the composition of our hospital changes over the course of a year.

plot_full_beds <- function(){
  options(gganimate.nframes = max(full_beds.df$Day)*2)
  full_beds.gg <-ggplot(full_beds.df) + 
                 geom_point(aes(2, -0.3,  size = Healthy), color = "royalblue") +
                 geom_point(aes(1,   1,   size = ED),  color = "deepskyblue") + 
                 geom_point(aes(2,   1,   size = ICU), color = "deepskyblue") +
                 geom_point(aes(3,   1,   size = Gen), color = "deepskyblue") +
                 geom_point(aes(2.5, 3,   size = Deceased), color = "steelblue") +
                 geom_point(aes(2,   2,   size = Transferred), color = "steelblue") +
                 geom_point(aes(1.5, 3,   size = Cured), color = "steelblue") +
                 
                 geom_text(aes(2,  -0.5, label = paste("Healthy\n", Healthy))) + 
                 geom_text(aes(1,   0.7, label = paste("ED\n", ED))) +
                 geom_text(aes(2,   0.7, label = paste("ICU\n", ICU))) +
                 geom_text(aes(3,   0.7, label = paste("Gen\n", Gen))) +
                 geom_text(aes(2.5, 2.7, label = paste("Deceased\n", Deceased))) +             
                 geom_text(aes(2,   1.7, label = paste("Transferred\n", Transferred))) +             
                 geom_text(aes(1.5, 2.7, label = paste("Cured\n", Cured))) + 
                 
                 transition_states(Day,
                                   transition_length = 2,
                                   state_length = 2) +
                 scale_size(range = c(1,80)) +
                 theme_bw() +
                 theme(axis.line=element_blank(),
                       axis.text.x=element_blank(),
                       axis.text.y=element_blank(),
                       axis.ticks=element_blank(),
                       axis.title.x=element_blank(),
                       axis.title.y=element_blank(),
                       legend.position="none",
                       panel.background=element_blank(),
                       panel.border=element_blank(),
                       panel.grid.minor=element_blank()) +
                 xlim(c( 0.5, 3.5)) +
                 ylim(c( -1,  3.5)) +
                 ggtitle("Day {closest_state}") +
                 theme()
  print(full_beds.gg)
  anim_save("Hospital_simulation.gif", full_beds.gg)
}

plot_full_beds()

Conclusions:

The base logic of the model has wide applications, and not just in the medical field. We could take the model far out of its original purpose, and with minor modifications this could, for example, serve as a way to sort and transfer children into different classes.

Overall, this model is a simple one, and that is intentional. Many possible additions come to mind when looking at this model. We could add doctors, nurses, and administrators. We could try different amounts of beds in each area. We could add budget constraints. We could have a waiting room, or have the simulation tick away at an hourly rather than daily basis, or try different distributions for severity change, or have the ability to set up temporary triage centers for emergencies, or… any one of a million different things.

But here, the biggest constraint is not things that could be added to the simulation. Instead, the problem is the usual culprit, not enough good data. Without accurate data, we would likely soon find that a more complicated simulation would make it impossible to form proper conclusions.